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SUMMARY 


Two matrix methods for calculating responses of linear systems for 
given forcing functions and indicial responses are analyzed in order to 
find a method suitable for solving the inverse problem of calculating 
forcing functions from known responses. The first method consists of 
a numerical evaluation of Duhamel’s integral, the second of an approxi- 
mation of the forcing function by means of straight— line or parabolic- 
arc segments. The second method is more generally adaptable to the 
solution of the inverse problem than the first; when the first method 
can be used, however, it is usually less time consuir' 

The results of both methods are compared with exact solutions for 
both the direct and the inverse problems. Both methods are found to be 
convenient to use and sufficiently accurate for many practical purposes. 
They may find application in work on gust loads, maneuvering loads, 
impact loads, in electric problems, as well as many other dynamic and 
some static problems. 


INTRODUCTION 


It is often of interest to calculate the response of a static or 
dynamic system to an arbitrary forcing function, such as the response 
of an airplane to an arbitrary gust or stick— force variation. Similarly, 
it is often of interest to calculate the forcing function that gives rise 
to a known response. 

If the system is linear so that solutions may be superposed, the 
direct problem of calculating the response to an arbitrary forcing 
function may be solved by superposition with the use of an indicial 
response (the response to a unit jump function, see fig. l) . The integral 
form of this superposition is known as Duhamel's integral. When the 
integral cannot be evaluated in closed form either numerical or graphical 
methods (reference l) may be used. 
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For the purpose of solving the Inverse problem of calculating 
forcing functions from known responses, Duhamel 8 s integral may be 
considered to be an integro-differential equation. Because solutions 
to this equation cannot always be obtained in closed form and graphical 
methods cannot conveniently be adapted to this problem, numerical methods 
must be resorted to. 

In this paper two numerical methods are presented for solving the 
inverse problem. One method consists of a numerical evaluation of Duhamel*s 
integral in a manner similar to that of Simpson^ rule. The other method 
consists of an approximation of the forcing function by means of straight- 
line or parabolic-arc segments. Both methods are expressed in matrix form 
and are baBed on matrices that are obtained by inverting matrices which 
constitute the solution to the direct problem. Consequently, the solution 
of the direct problem by the two matrix methods is discussed first and the 
adaptation of these methods to the solution of the inverse problem, later. 
The mechanics of computing solutions to specific problems are discussed in 
some detail in the corresponding sections of this paper concerned with the 
analysis of the direct and inverse problems. Methods of evaluating the 
responses to unit gradient and unit parabolic forcing functions (see fig. l) 
required in the second method of this paper are described in the appendix. 

The adaptation of these and similar numerical methods to the solution 
of the other inverse problem, that of calculating Indicial responses from 
known forcing functions and responses, is discussed briefly. 


SYMBOLS 

s independent variable, usually time 

A indicial response 

G arbitrary forcing function 

B response to arbitrary forcing function 

B response to unit gradient function 

C response to unit parabolic function 

f integrand of Duhamel , s integral 

cr variable of integration corresponding to s 

| independent variables equivalent to s 
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differentiating matrix 

matrix which, performs the function of Duhamel’s integral 

constants of linear and parabolic term in the parabolic— arc 
appr oximat i on 

nth point on any curve or number of last point of interest 

coefficient in Q or Q^_ matrix or coefficient of original 
matrix 

b coefficient in inverse of Q or Q]_ matrix or coefficient 

in auxiliary matrix 

c,d,e,f coefficients in Qg matrix 
g submatrix of matrix 

h submatrix of inverse of Q,g matrix 

Matrix notation: 


M 

M 

n 

a 


[] 

{} 


square matrix 
column matrix 


The subscripts on the k, B, C, G, and R functions or their 
derivatives denote the point at which the function (or derivative) 
is taken; thus, for instance, G* = G ! (s) . 


A prime mark indicates differentiation with respect to s or a. 
An underscore indicates a submatrix. 


ANAITSIS OF THE DIRECT EROBLEM 
Numerical Evaluation of Duhamel s s Integral 


The response R(s) of a linear static or dynamic system to an 
arbitrary forcing function G(s) can be obtained from the in dicia! 
response of the system A(s) by superposition. If the forc ing 
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function is approximated "by a staircase function (fig. 2 ) } the response 
at different values of the independent variable s will be: 


K(0) 

ECsj) 

E(s 2 ) 

K(s 3 ) 


E 0 ~ A 0 G 0 


E i ~ A i^o + ^0(^*1 

e 2 ~ 

E^ £ A^Gq + a 2 (g*j^ 


G-q) 

Gq) + Aq(G2 
Gq) + A^(Gr> 


- Gj_) 

— G-^) + Aq(G^ — G 2 ) 


> 


( 1 ) 


As the interval As approaches zero these discrete values of E 
approach a function E(s) given by Duhamel*s integral: 


E(s) = A(s)Gq + A(s - cr)G*(cr)dcr (2) 


dO 


where a is the variable of integration corresponding to s. 

For cases in which Gq is not zero the integral in equation (2) 

can be evaluated as indicated in the following analysis | to this integral 
the function A(s)Gq may then be added. When Gq is zero 5 equation (2) 
reduces to the simpler form: 


E(s) = A(s — cr)G* (cr)dcr (2a) 

This expression is evaluated numerically in two steps: an expression 

for the derivative of the forcing function G(a) is first determined 
and then the integration of the function A(s — c)G 3 (a) is performed 
by Simpson* s rule. 

An approx imat e value of the derivative of G at point Sn may be 
obtained from the relation 


QJ <v 

^ n 


G n + l 




As 


( 3 ) 
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which, in effect implies that the actual curve for G is replaced by a 
parabolic segment between the points s n _l and s n| l. If, in addition, 

the slope at s = 0 is estimated as 


G *0 


G l/2 

As/2 


(3a) 


a differentiating matrix for calculating G*(s) can be set up from 
equations ( 3 ) and (3a): 


p’o 


~~ 2 0 0 0..' 


°T /2 

G ’l 


-1 1 00... 


G3/2 

< G*2 

- 

f ^ A 

As 

0-1 10... 



G *3 


0 0-11... 


Grj/2 

- V 




• 


(4 ~ ST® { G } 


The integration of the function A(s — cr)G* (ct) = f(cr) indicated in 
equation (2a) may be performed by Simpson* s rule. For even values of n, 

P°h 

R n = J f(o)da 

~ Aj^(fQ + 4f p + 2f2 + 4f^ + 2fj^ + . . . + ^f n _p + f n ) 

For odd intervals the integration over the last segment s- n _-[_ to 
can be performed by passing a parabola through the three points s n _2^ 
s n _ 1 , and s n and integrating the parabola between s^p and s n . 

This procedure yields a set of coefficients (— l/l2, 2 / 3 , and 5/l2) which 
is used in the same manner as the Simpson integrating factors ( 1 / 3 , 4/3, 
and 1 / 3 ) which would be obtained upon integrating from s n _2 to s n . 


(5) 

8 n 
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Thus for odd values of n, 


Rn = 


O^n 


f (cr)dcr 


UO 


~ ^(fo + hfq + 2f2 + 4 f 3 + 2f k + • 

+ As^f n _ 2 + f n _ x + ^f n ^ 


+ 4fn-J+ + 2fjx— 3) 


(5a) 


Since there are only two values of the function f that can he used to 
calculate Rq, the trapezoidal rule is used in the initial interval 
instead of Simp 8 on*s rule, so that 


% ~ f( f i + f o) 


( 5 "b) 


Equations (5) , ( 5 a) , and ( 5 h) can he combined hy writing them in 
matrix form: 


or 


% 


Eg 


< E 3 

As 



R 5 ! 






§Ao 

0 

0 

f 2 

b 

3^° 

0 

— Ao 

3 3 


Al 



h 

|*2 

b 


b 

^*3 
3 3 

h 



{*} 

X As Jj 


0 

0 

0 


0 


0 

0 

0 

0 

12 


0 


G' 


G* 


G 1 


0 


L. 


( 6 ) 


[A}{ G ^ (6a) 

Then, substituting the expression for -|g s j" from equation (ha) into 
equation (6a) gives 

{r} ^ [ A] |d]{g} (6h) 

= [Q]{g} ( 6 c) 

where jo] is a combined matrix that performs the operations indicated by 
Duhamel 1 s integral numerically. 
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The steps required in computing responses by this method can he 
summarized as follows: 

1. The values of the forcing function are tabulated for the midpoint 
of each interval and those of the indicial function are tabulated for 
each end point. 

2. The values of the forcing function are then multiplied by the 

D matrix (equations A) and (^a) ) to obtain the values of a matrix that 
consists of As times the values of G* (the derivative of the forcing 
function). If n is the last point of interest, the tabulation of the 
values of the forcing function should be extended to the value at the 
midpoint of the next interval (3^ ,1: the values of G* can then be 

calculated up to G* n . 

3 . The values of the indicial response are multiplied by the 
factors l/2, l/ 3 , 5/12, 5A, and so forth and tabulated in the manner 
indicated by the A matrix of equation ( 6 ). The A matrix has one 
more column than it has rows. If n is the last point of interest, 
it will have n + 1 columns, which correspond to the n + 1 values 
of G 1 (starting with G' 0 ). 

4. The A matrix is postmultiplied by the column As {g^ obtained 
previously. The values calculated in this manner are the values of the 
response at the end points of the intervals, provided Go is zero. 

Only n values of R are calculated; Rq is zero. 

5. If Gq is not zero, the response is obtained by adding the 
terms GqA-^, GqA 2 , GqA-^, • • • to the values of Rj_, Rg, R 3 , . . . 

calculated in the preceding step. The initial value of the response 
is then GqAq . 


Approximation of the Forcing Function 

In the preceding section the response to an arbitrary forcing 
function was evaluated numerically by approximating the integrand 
A(s — ct)g 8 (ct) of Duhamel*s integral by parabolic segments . Another 
method of calculating such a response consists of approximating the 
forcing function either by straight-line segments or by parabolic-arc 
segments. The response is then calculated as the linear superposition 
of responses to the unit gradient and unit parabolic forcing functions 
shown in figure 1 . 

Straight-line appro ximation.— If the arbitrary forcing function G(s) 
is approximated by straight-line segments as shown in figure 3 , the values 
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of the response to that forcing function at different values of s c an 
he determined from the response B(s) to a unit gradient function (see 
fig. 4) as follows; 


*0 = 0 


E 1 ~ 

Ep ^ Bp G-^ + ( Gp — G^_ ) 

R3 ~ B3G1 + Bp(Gp - G^) + B-j^G-^ — G2 ) 




J 


This relation can he expressed in matrix form as 


( 7 ) 


or 


f- 

% 


Ep 


I * 3 

> ~ 

E 4 

e 5 





B x 0 0 0 0 

Br) 0 0 0 

B ^ B/p B^ 0 0 

Bj^ B^ Bg Bj 0 
B^ B]^ B^ B-p Bj 


{e} 3 [b][di]{g} 
- 


1 0 0 0 0 ... 


°1 

-1 1 0 0 0 ... 


Gp 

0-1 1 0 0 ... 

< 

g 3 > 

0 0-1 10 ... 



0 0 0 -1 1 ... 



« <* a t • • 0 • 




(7a) 


(7b) 

(7c) 


A method for calculating the response to the unit gradient function 
B(s) in terms of the indicial response A(s) is presented in the appendix. 
The steps required in computing responses by this method can then he 
summarized as follows: 


1. The values of the function G- are tabulated at the end points 

of the intervals from point 1 to point n 3 and those of the function k, 
at the midpoints as well as at the end points from point 0 to point n. 

2. The values of the function B at the end points of the intervals 
(from points 1 to n) are calculated' from equation (A 6) of the appendix. 
The matrix of equation (A 6) will have n rows and 2n + 1 col umn s, 
since 2n + 1 values of A are used (at the n end points., at the n 
midpoints * and at 0) and n values of B are calculated. 
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3 . The values of the function B calculated in this manner are 
entered in a matrix as shown in equation ( 7 a) and are postmultiplied 
"by the Dp matrix to yield the Qp matrix. Both the Dp and the Qp 
matrices are square and of order n. 

The response to the forcing" function G(s) is found by post- 
multiplying the Qp matrix by the G matrix as indicated in equation (7c), 
provided Gq is zero. 

5. If Gq is not zero, it must be subtracted from every other value 
of G in tabulating the G matrix. The terms GQAp, GqAq, . . . must 
then be added to the values of the response calculated in the preceding 
step for points 1, 2, ... . The response at point 0 is then GqAq. 

Parabolic-arc approximation .— A closer approximation to the forcing 
function than the straight-line approximation of the preceding section 
may be had by fitting parabolic-arc segments to the function. In the 
first segment 0 < s < As, for instance, the forcing function may be 
expressed approximately as 


G(s) 2M -if-*-) 2 + N n 

\&3J x \Asy 


( 8 ) 


provided Gq is zero. The constants Mp and Np are selected in such 
a manner that the parabolic arc coincides with the true values of G at 

s = 0, s = — , and s = As; therefore. 


Mp = — hGp Iq + 2Gp 

*1 = kG l/2 ~ Gl j 

similarly, in the second interval 

G(s) - °i~ + 


where 


M 2 = -KG 3/2 ~ Gp) + 2(02 ~ Gp) 

N 2 = K G 3/2 ~ ^ _ (p2 ~ Gp) 


( 9 ) 


(8a) 


(9a) 
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and in further intervals, for s n _i = s = Bq, 

G(b) - ^-1 * + 

-where 

•t - -*►(<*- 1 - %- l ) + 2 ( B tt - Oh - l ) 
Mn = 4(s_l - %_i) - (Oh - Uhl) 

These relations may he written in matrix form as 


m-l 


-4 

2 

0 

0 

0 

0 . 



Gl/2 

% 


k 

-1 

0 

0 

0 

0 . 



°1 



0 

2 

-4 

2 

0 

0 . 

9 9 


G 3 / 2 

n 2 

> = 

0 

-3 

k 

-1 

0 

0 . 

9 9 

■4 

Gip ** 

m 3 


0 

0 

0 

2 -4 

2 . 

9 *9 


°5/2 

n 3 


0 

0 

0 

-3 

h 

-1 . 

9 9 


°3 

• 


— 

* 

• 

» 

* 

• 9 

9 9 


9 


(8b) 


(9b) 


( 9 c) 


or 

{m} = |dJ] {g} ( 9 d) 


Once the M and N coefficients pertaining to the arbitrary 
forcing function G(s) have been evaluated, the response to that forcing 
function at different values of s may be determined from the 
response B(s) to a unit gradient function and the response C(s) to 
a unit parabolic function as follows: 
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or 


or 


*1/2 ~ C l /^1 + B l/2%. 

^ + B^N*^ 

^3/2 ^ ^3/2^1 ^ ■^3/2^'X ^1/2^2 ^ y/2^2 

Eg /v CgM-L + BgN-L + C-jMg + B-jNg 

^ 5/2 ^ ^ 5 / 2^1 ^ 5 / 2^1 ^ 3 / 2^2 ^ ^ 3 / 2^2 ^ ^ l / 2^3 ~^~ 


E l/2 


B^ ^g 0 0 0 0 ... 


% 

% 


C 2 B x 0 0 0 0 ... 


N 1 

E 3 /2 


C3/2 B 3 /g C-L/g B-L/g 0 0 . . . 

-c 

M2 k 

Eg 


C2 Bq c x b x 0 0 • 1 • 

1 

n 2 

* 5/2 


C 5/2 B 5/2 C 3/2 B 3/2 C l/2 B l/2 * * * 


m 3 



• • 0 • • 1 t * • 


* 


{e> X [c] {m} 




(10) 




(lOa) 


(lOb) 


where the M^N column ia calculated by means of equation (9c) so that 

( E ) ~ [°] M { G } (10o) 


which may also be written as 



(lOd) 


where the Qg matrix in effect performs the operations of Buhamel 8 s 
integral. v ... 
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If values of R are desired at only integral multiples of As, 
alternate rows in equation (10a) may be omitted so that the equation 
is simplified to 


% 


C-l B x 0 0 ' 0 0 ... 


*L 

>2 

> ~ 
A# 

C 2 B 2 Cl Bq 0 0 ... 

< 

Nil 

r 3 


B^ Cg Bg B^ . . . 


Mg 

« 


Lr Jj 


Isl 


(lOe) 


Methods for calculating the responses B(s) and C(s) in terms of 
the indicial response A(b) are presented in the appendix. The steps 
required in computing responses by the method of parabolic-arc approxi- 
mation to the forcing function may then be summarized as follows: 

1. The values of the forcing function as well as those of the 
indicial response A(s) are tabulated at the end points and midpoints 
of the intervals] there will then be 2n values of G (from Gj/g 

to G^) as well as 2n + 1 values of A (from Aq to A^) . 

2. If values of the response are desired only at the end points of 
the intervals, the functions B and C can be calculated at those 
points from equations (a 6) and (Al^ ) (disregarding alternate rows). If, 
on the other hand, values of the response are also desired at the mid- 
points, the functions B and C must be calculated at those points as 
well; equations (A7) and (Alh) may be used for this purpose. Equations (A8) 
and (A15) may be used to yield somewhat more accurate values of B 

and Cq /2 than furnished by equations (Aj) and (Al^), respectively. 

3 . The values of B and C are tabulated in matrices as indicated 
in equation (lOe) or equation (10a), depending on whether the response is 
to be calculated only at the end points or at both the end points and the 
midpoints. 

4. The values of the coefficients M and N are calculated from 
the values of G by means of equation (9) . 

5. The response to the function G is calculated from equations (10a) 
or (lOe), provided Go is zero. If Go is not zero, a term 

(GqA^, GqA 2 , . . .) has to be added to each value of the response 

calculated from equation (10a) or (lOe); the response at s = 0 is 
then GqAq . 
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ANALYSIS OF THE INVERSE PROBLEMS 


There are two types of inverse problems. The one concerned with 
calculating the forcing function from a known response and indicial 
response appears to he generally of ftore interest. It is consequently 
treated in more detail in this paper than the second inverse problem, 
which consists of calculating indicial responses from known forcing 
functions and responses. 


Calculation of the Forcing Function 

For the purpose of calculating responses to arbitrary forcing 
functions, the numerical methods presented (equations (6c), (7c), 
and (lOd)) may be considered to consist of summation formulas involving 
the indicial response, the arbitrary forcing function, as well as 
certain numerical factors. For the purpose of solving the inverse 
problem, however, the same equations may be considered to be linear 
simultaneous equations with known responses on one side and unknown 
values of the forcing function multiplied by certain coefficients on 
the other side. These coefficients are the elements of the Q matrices. 
The unknown values of the forcing function may then be obtained by 
solving the simultaneous equations in any convenient manner. There must, 
of course, be as many equations as unknowns, or, in other words, the Q 
matrices must be square if they are to be used for the solution of the 
inverse problem. 

In the following analysis Gq is assumed to be zero. When it is 
not zero it may be calculated from the relation 


G -5fi 

G °'a 0 


(ID 


and the terms GqA-^, GqAp, . . . must be subtracted from the known values 

of the response before they are operated upon in the manner indicated in 
the following sections in order to calculate the forcing function. To 
the values of the forcing function calculated in this manner the value Gq 
as calculated from equation (ll) must then be added at each point. 

Numerical— integration method .— The numerical evaluation of D uhame i 9 s 
integral leads to equation (6c). In the case of the inverse problem, a 
unique solution for the values of G corresponding to a given set of 
values of R can be obtained from this equation only if A q is zero, 
since otherwise there are more unknown values of G than there are 
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equations. Consequently this method can he applied to the analysis of 
the inverse problem only if A q is zero. In that special case the Q 
matrix is triangular j the solution of equation (6c) can be performed by 
the following method: The Q matrix may be written in the form: 


'Q 


a ll 

0 

0 

0 

a2i 

CVJ 

0 ^ 

0 

0 

a 3 l 

a 32 

a 33 

0 

& kl 

a k2 

a 4 3 

a A4 


( 12 ) 


where the coefficients a are obtained by postmultiplying the A 
matrix of equation (6) by the D matrix of equation (k) . The solution 
of equation (6c) is then 


°l/2 = iJI E 1 

G 3/2 = S^( R 2 ~ a 2l9l/2) 

<b/2 = a^( E 3 “ a 31 G l/2 “ a 32 G 3/2) 


If many response functions are to "be analyzed for the same indicial 
response s it may he expedient to invert the Q, matrix as follows! 



t-UL 0 0 0 

h^i ^22 0 0 

"^31 ^32 ^33 ® 

bin Hz ^3 


( 1*0 
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where each column of jo] - 1 is obtained successively by means of 
equation ( 13 ) from an appropriate R column of the type 



and the values of G pertinent to the given problem. 

The steps required in calculating forcing functions by this method 

are: 

1. The n values of the response are tabulated at the end points 
of the intervals from points 1 to n. 

2. The A matrix is obtained from equation (6) as described in the 
section concerned with the direct problem; it is postmultiplied by the 

D matrix of equation (4) to yield the Q matrix. 

3 . The elements of the Q matrix are considered to be the coeffi- 
cients of simultaneous equations for the unknown values of G from G^y, 

to G n _i in terms of the known values from E-j_ to R n . The solution 
2 

of these equations may be carried out as outlined in equation ( 13 ) . 

k. If it is desired to invert the Q matrix, equation (l4) can be 
used. Tne values of the forcing function corresponding to any set of 
values of the response may then be obtained by premultiplying the values 
of the response by the inverse of the Q matrix. 

Straight-line approximation .— The straight-line approximation to 
the forcing function leads to equation (7c), which involves the 
triangular matrix 


a x 0 0 

So a-j_ 0 


[ 8 .^ 


0 

0 

0 



a^ a^ ag a^_ . . . 


( 15 ) 
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■where 

a l “ B 1 

= B 2 

a 3 = b 3 



When there are n values of G (from Gq to Gn) and n values of R 
(from Rq to R ^) } the Qq matrix is square and there are n equa- 
tions for the n unknown values of G. 

The solution for these values is carried out very easily as a 
result of the triangular form of the Qq matrix. In fact 


°1 " a^ E l 

“ ^( E 2 “ 


g 3 

% 



a^Gq — agGr^ 


— aj^Gq — a^Gp 



( 17 ) 


If the forcing functions for many known responses are to he obtained for 
the s ame system (with the same indicial response it may be more expe- 
dient to calculate the inverse of the Qq matrix than to solve for each 
forcing function separately. The inverse of the Qq matrix may be 
written as 


bp 0 0 0 

b 2 bp 0 0 

b^ b 2 bp 0 
blp *>3 b 2 bp 

L* 


-1 


(18) 
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where 



3 = - a7( a 3 b l ~ a 2^2^ 


"bo = 


h^ = - ^(aj^bq - a 3 h 2 - a 2 h 3 ) 



(19) 


The steps required in calculating forcing functions by this method 

are: 

1. The values of the response are tabulated at the end points of 
the intervals from points 1 to n. 

2. The Qq matrix is calculated as indicated for the direct 
problem. 


3. A set of simultaneous equations having the elements of the Qq 
matrix as coefficients and the values of R as knowns is solved for 
the unknown values of the forcing function G]_ to as indicated 

in equation (17). 

If it is desired to invert the Qq matrix, equations (l8) 
and (19) can be used; the forcing functions may then be obtained by 
premultiplying the given sets of values of R by the inverse of the 
Qq matrix. 

Parabolic-arc approximation .— The Q 2 matrix obtained for the 
parabolic-arc method (equation (lOd)) is not quite triangular but may 
be reduced to triangular form by partitioning: 



18 


RAC A TR 1965 


dLj | 0 0 j 0 0 ■ . . 

©1 "^1 | ® 0 j 0 0 • e a 



C£ &2 I c ± d.^ loo • * • 

©2 f*2 j ®1 -£*1 1 0 0 • . • 

| - 
°3 d 3 j c 2 *2 °1 d l • • • 

©3 f 3 I ©2 f 2 e l * • • 




g>2 ^ « a * 


83 §2 Si • • • 


( 20 ) 


where 



ci d x 

f 1 f l 


82 


C2 d.2 
e 2 f 2 




(20a) 



0 

0 


and where the c, d, e, and f values are obtained by postmnltiplying 
the C matrix of equation (10b) by the Dg matrix of equation (9d) to 
yield the Qg matrix. If the column matrix of response values is also 
partitioned into pairs such that 
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E l/2 


E n 


{*} =Je 3/2 ] 
*2 


% 


( 21 ) 


•where 


R 1 


E 2 




V 

S 

=1 ^ 


E 3 / 2 


*2 

A 




(21a) 


then a solution of the simultaneous equations having the elements of 
the Q 2 matrix as coefficients may he effected in the manner indicated 
for the straight-line method, except that two— by— two matrices take the 
place of ordinary algebraic quantities. Trie g matrices of 
equations (20a) take the place of the coefficients a (equations ( 15 ) 
and (l6)), the E matrices of equations (21a) take the place of the 
values of E in equation (17), and the forcing function is computed in 
the form of G matrices defined by 
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which take the place of the values of G in equation (17). If the 
Q2 matrix is to he inverted, the resulting matrix has the form 


h-L 0 0 ... 



hg hj_ 0 ... 

hg h 2 hj . . . 



(23) 


where the h matrices take the place of values of t in equation (18) 
and are calculated as indicated in equation (19). The quantity l/aq 

T 

is replaced hy gq , the inverse of the gq matrix. The operations 
called for in equations (17) and (19) are now, of course, matrix multi- 
plications, so that the order of the multiplicands must he preserved. 
The fact that the inversion of the Qg matrix entails, with the 
exception of the inversion of the two— hy— two gq matrix, only matrix 
multiplications and additions of low— order matrices facilitates the 
problem greatly. 

The steps required to calculated forcing functions hy this method 

are: 


1. The values of the response are tabulated at the end points and 
midpoints of the intervals from Rq/2 to Rq. 

2. The values of B and C are calculated at both the end points 
and the midpoints of the intervals and entered in a matrix as shown in 
equation (10a) . This matrix is postmult iplied hy the Dg matrix shown 
in equation (9c) to yield the Q2 matrix. 

3. A set of simultaneous equations having the elements of the 

Qg matrix as coefficients and the values of R as knowns is solved 
for the unknown, values of the forcing function to in the 

manner suggested in equation (17), where the submatrices indicated in 
equations (20a), (21a), and (22) are used instead of the terms a, R, 
and G of equation (17). 

4. If it is preferred to solve the problem hy inverting the Q2 
matrix, the inversion may he performed as indicated in equations (l8) 
and (19), where the submatrices indicated in equation (20a) are used 
instead of the values of a of equation (19), however) the forcing 
functions may then he obtained hy premultiplying the given sets of 
values of R hy the inverse of the Q2 matrix. 
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Calculation of the Indicial Response 

The numerical— integration method outlined for the direct problem 
may in certain cases he modified to obtain a solution to the second 
inverse problem of calculating an Indicial response from a known forcing 
function and response. Equation (6) may be rewritten in the form: 



which permits a solution for the values of A from known values of R 
and G® . The values of G’ may be. obtained from the values of G at 
the midpoint of the chosen intervals by means of equation (4) . As in 
equation (6), it is assumed that Gq is zero. Furthermore unless 
either Aq or Go is zero as well, equation (24) will have more 
unknown values of A than equations, so that it cannot be solved. 

If these conditions are not met, it may be possible to use another 
approach which consists of using the alternate form of Duhamel s s integral: 


R(s) = A 1 (cr)G(s — a)dcr (25) 

l 0 


The following equations may be obtained in a manner analogous to that 
used in obtaining equation (6) : 




(26) 
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This equation serves to solve for the values of A® in terms of the 
known values of the response; the values of the forcing function appear 
as coefficients in the simultaneous equations, as do the integrating 
factors. Equation (25) is valid when A q is zero; it does not imply 
any restriction on Gq. However, unless either Go or Ao is zero 
there will he more unknown values of A than equations in equation (2 6), 
so that the equation cannot he solved. A set of values of A may he 
obtained from the calculated values of A* hy means of an integrating 
matrix; 


A £ 

A2 


a 3 

> = As 

A k 

b 

• 



5/12 

1/3 

5/12 

1/3 

5/12 


2/3 

ty3 

1 

V3 

1 


-1/12 

0 

0 

0 

1/3 

0 

0 

0 

5 A 

1/3 

0 

0 

2/3 

V3 

2/3 

0 

5 A 
• 

2/3 

e 

V3 

• 

i/3 

0 



( \ 


A*0 


A ’l 

< 

A *2 > 
|A* 


3 


\ a \ 


l • > 


(27) 


The integrating process indicated in equation (27) will tend to 
average out any discrepancies in the values of A* calculated from 
equation (2 6) . 

Since in the particular cases where equations (2k) or (2 6) may he 
used the matrices of the values of G or G® will be triangular, the 
methods indicated in the preceding sections in analyzing the numerical 
integration method will also he applicable to the solution of these 
equations or the inversion of these matrices. 


RESULTS AMD DISCUSSION 
Comparison of Numerical and Exact Results 


The direct problem .— The accuracy of the numerical results has been 
investigated for several cases hy comparing their results with known 
exact solutions for simple indicial response and forcing functions. 

The indicial response selected is; 


A(s) = 1 


1 —0.2s 
— e 

2 


(28) 
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Two forcing f unctions have "been considered 


and 


Gvj-(s) = sin -- B 


( 29 ) 


(TnCs) 



COS 



( 30 ) 


The responses ho the two forcing functions have been computed by 
Duhamel’s integral exactly and by the numerical methods of this paper. 

The interval was chosen as As = 2, which is one twenty— fourth of the 
period of the forcing function. The response B(s) to the uni t gradient 
and the response C(s) to the unit parabolic function required in the 
second method (the straight-line or parabolic— arc approximations) have 
been computed by the numerical methods indicated by equations (A7) 
and (Aik) ; they are presented in table I with the exact values. , A 
comparison shows that the numerical methods for the evaluation of the 
unit responses B(s) and C(s) should yield results to an accuracy of 
better than 0.1 percent for reasonably small intervals . 

The response to the forcing function of equation (29) has been 
calculated by all three numerical methods and is presented in table II. 

A comparison of the results indicates that the numerical— integration 
method and the straight-line approximation have the same accuracy, 
approximately 0.2 percent. The parabolic-arc approximation gives much 
better accuracy, approximately 0.01 percent in this example. 

The inverse problem .— The forcing functions which give rise to the 
responses may be calculated from these exact responses and the indicial 
response given by equation (29) - Since A q is not zero, the approximate- 
integration method cannot be used. Therefore, only the other two methods 
have been used to calculate the forcing functions. 

An increment As = 2 was used and the required unit gradient and 
unit parabolic responses were taken from table I. The forcing functions 
calculated by these methods should check those given exactly by 
equations ( 30) and ( 31 ). The comparison is presented in table III. 

The average accuracy of the straight-line approximation applied to the 
inverse problem is approximately 0.3 percent; that of the parabolic-arc 
approximation, approximately O.Oh percent. No points have been calcu- 
lated at the midpoints of the intervals by means of the straight-line 



2k 


NACA TN 1965 


approximation since they would simply he the average of the ordinates 
at the end points of the given intervals in view of the approximation 
inherent in that method. 


Factors Affecting the Accuracy of the Methods 

The comparisons of the numerical methods with known solutions 
show that, with intervals of the order of l/20 to l/30 of the period 
of the first natural frequency of the system or that of the forcing 
function, the numerical— integration method and the straight-line 
approximation yield results that have an accuracy of 1 percent or 
hetter; whereas the results of the parabolic— arc approximation are 
accurate to 0.1 percent or better. Similar calculations have indicated 
that an accuracy comparable to these values may be expected even when 
higher harmonics are present, provided those higher than the third are 
unimportant compared with the first few. 

In general, the accuracy of the numerical methods depends to some 
extent on the shape of the curve that is approximated; the closer the 
approximate curve (which consists of straight-line or parabolic-arc 
segments) fits the actual curve, the higher the accuracy. Consequently, 
the interval should be chosen smal 1 enough that the degree of curvature 
of the approximated curve is small within a segment; any reflexes in 
the approximated curve should be near end points of segments, if 
possible. 

In the case of the inverse problem, great care must be exercised 
in the case where both the indicial response and its first derivative 
vanish initially (Aq = A*q = 0). The determinant of the coefficients 
of the simultaneous equations to be solved for the values of the 
forcing function tends to be relatively sma ll for this case and the 
equations relatively ill— behaved. It is then difficult to obtain 
reliable results by numerical means unless very small increments are 
taken near s = 0. The accuracy of the values of the forcing function 
may, however, be Improved by replacing the actual indicial response by 
a straight-line segment extending over the first two or three intervals 
for the purpose of computation; in this manner a nonvanishing value is 
assigned to A*q and the simultaneous equations for the values of Gl- 
ean be solved more accurately. 


Factors Affecting Choice of Methods 

The accuracy of the numerical— integration method and straight-line 
approximation for both problems tends to be approximately the same. The 
numerical— integration method is less time consuming than the other method 
if only one response or one forcing function is to be analyzed, since the 
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straight-line approximation requires the computation of the function B. 
If many cases are to be analyzed for the same indicial response, the 
straight-line approximation will tend to "be the more expedient, since 
it involves matrices which are somewhat more convenient to use. The 
numerical-approximation method is, of course, applicable to the inverse 
problem only when Aq is zero. 

The parabolic-arc approximation is more accurate than the other 
methods but is also more time consuming. If by decreasing the size 
of the segments used in the other methods their accuracy is increased 
until it is comparable to that of the parabolic-arc approximation, the 
expenditure of time required becomes comparable as well. The shape of 
the function to be approximated will then determine whether any small 
advantage would be with the parabolic-arc approximation or the other 
methods . 


CONCLUDING REMAKES 


Two matrix methods for solving the problem of calculating forcing 
functions from known responses have been derived by means of a numerical 
analysis of the problem of calculating responses fro.,. -cing 

functions. The first method, which consists of an approximate evalu- 
ation of Duhamel ? s integral, leads to a matrix that can only be inverted 
when the initial value of the indicial response is zero. The matrices 
obtained in the other method, which consists of approximating the 
forcing function by straight-line or parabolic-arc segments, can always 
be inverted readily. Some results obtained by the numerical methods for 
simple explicit functions have been compared with exact results. The 
accuracy of the numerical methods was found to be adequate for most 
practical purposes. A brief discussion has been given of the possi- 
bility of calculating the indicial response by the approximate- 
integration method If the forcing function and the response to it 
are known. 


Langley Aeronautical Laboratory 

National Advisory Committee for Aeronautics 

Langley Air Force Base, Ya., July iL, 19^9 



2 6 


NACA TN 1965 


APPENDIX 

EVALUATION OF EESPONSES TO UNIT GRADIENT 
AND UNIT PAEABOLIC FOECING FUNCTIONS 

The response to the unit gradient function can he determined from 
the indicial response by means of Duhamel’s integral. For the unit 
gradient , for 0 < s % As, 


Gr(s) 


a 

As 


and. 


G*(s) = 


As 


for s As, 


and 


G(s) = 1 


G* (s) = 0 


so that, for 0 ^ s ^ As, 


B(s) = 

As 


fls 


do 


A(s — cx)da 


(Al) 


and, for s > As, 


B(s) 


_1_ 

As 


PAs 


A(s — cx)dcr 


do 


(A 2 ) 
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or, if £ is substituted for a — a and I is substituted 
for — (b — As), for As* 


B(s) 



P S A(?)d? 


10 


(A3) 


and, for s n > As, 



•As 


A( s n~ 


l)dl 


(a 4 ) 


If equations (A 3 ) and (a 4 ) cannot be evaluated in closed form, 
numerical methods can be used. Since the accuracy of the further calcu- 
lations (equation ( 7 a)) depends on the accuracy with which the function B 
has been calculated, increments As /2 should be used in performing the 
integration indicated by equations (A 3 ) and (a 4 ). However, the values 
of B used need only be calculated at integral mull le r As. If 
Simpson’s integration rule is used to evaluate equation (A 4 ), the 
following values are obtained for B: 


or 


Bl ~ 3 (^0 + ^ A l/2 + A l) 

B2 ~^l( A l + 4A 3/2 + A 2) 


-"N 

B 1 


14 1 0 0 0 0 


f 

h 

^ A/ 1 

0 0 14 10 0 


A l/2 

b 3 

a 

6 

0 0 0 0 1 4 1 


A 1 f 

a 3/2 

« 


(A5) 


(A6) 


If for some reason the values of the function B are needed at the 
midpoints of the intervals, they can be calculated in a similar manner. 
The value at As /2 poses a special problem in that the Simpson factors 
are not directly applicable. However, the integration of function A 
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between 0 and As/2 may be carried out approximately by replacing A 
by a parabolic segment between 0 and As but integrating the parabola 
only between 0 and As/2, This procedure leads to the factors 5/12, 
2/3, and — l/l2 which are used instead of the Simpson factors l/3, *4-/3} 
and l/3 • Thus, 


B l/2 

B 1 

J B 3/2 

B 2 

b 5/2 


5/4 2 -l/k 0 0 0 



1 *4- 10 0 0 


A l/2 

0 1 A 10 0 


A 1 

00 1 A 1 0 

J. 

a r 
a 3/2 

0 0 0 1*41 



* 0 O ft ® ® 


L ■> 


(at) 


Should a more accurate -value of 
obtained from 


B l/2 


be desired* it can be 


B l/2 ~ |(A0 + *4Al/A + A 1 / 2 J 


r\j 

A; 


^(Aq + ^ A l/4 + A l/2j 


(AS) 


The response to the forcing function as approximated by equation (8) 
is composed of two parts. One part is due to a unit gradient B(s) which 
has been described in the preceding paragraphs. The second part of the 
response, that due to the unit parabolic function C(s), may similarly 
be calculated from the indicia!, response. For the parabolic function, 
for 0 ^ s < As, 


and 


for s > As. 
— 2 



G* (a) 


2s 

(As) 2 


G(s) = 1 
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and 


G*(s) = 0 


so that, for 0 ^ s < As, 


C(s) 


2 P B 

5 - A(s - 0)0 da 

(As) d do 


(A9) 


for s > As , 


0(b) 


df 


lAs 

I 

UO 


A(b — a) a da 


(A10) 


or, if £ is substituted for s-a and I is substituted for -(a -As), 
for 0 ^ s ^ As, 



If equations (All) and (A12) cannot be evaluated in closed form, 

numerical methods may be employed. The factor 1 — takes the 

As 

values 1 , l/ 2 , and 0 vhen the term s n - a is Sg^, s n _l, and s n , 

respectively, as does the factor — — & — for a = As. However 

As As * 

for s = ^ the factor -f- assumes the values l/ 2 , 0 , and - 1/2 

vhen a is As/2, 0, and -As/2; the point at a = is fictitious 

since it lies beyond the limits of the integration and is used only to 
furnish a better description of the curve in the region of interest. 
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These values may he combined -with the integrating factors to yield 


2_ 


J l/2 ~ As 2 


'1 ~ 


As 2 


(S ¥° + I (0)a i/ 2 - ik(- i) A i) 

^1/2 + (°) a 3^ 

+ i^a-l + (o)a 3 / 2 ^ 

°2 ~ As ^ 3 ^! + ^3/2 + 


C 0 / 

As 2 


^(^1/2 + ^ A 1 + 


> 


or, in matrix form. 


C l/2 


5/801/800. . . 


AO ^ 

C 1 


1 2 0 00... 


A l/2 

^ c 3/2 

r _ 

0 1 2 00... 

-c 

A 1 . 

C 2 

3 

• 

« 

• 

0 

OJ 

H 

O 

O 


a 3/2 

c 5/2 


0 0 0 12... 


Ag 

4 


• • « * • • * 




(A13) 


(A14) 


Again, as in equation (a 6), the rows in equation (Al4) for the values 
of C at the midpoints of the intervals may he disregarded if these 
values are not of interest. 

A slightly more accurate value of than that given in 

equations (A13) and (Al4) is given hy 

C l/2 ~ ^ if ^0 + ^ A lA + (° )A l/2) 


rv 

/v 


+ ^lA) 


(A15) 
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TABLE I.- COMPARISON OF THE NUMERICAL AND EXACT 
CALCULATIONS OF THE RESPONSES TO TEE UNIT 
GRADIENT AND UNIT PARABOLIC FUNCTIONS 


s 

Response to unit 
gradient function^ 
B(s) 

Response to unit 
parabolic function, 
C(s) 


Exact 

Numerical 

Exact 

Numerical 

0 

0 

0 

0 

0 

1 

• 5462 

.5468 

• 5317 

• 5317 

2 

.6290 

.6290 

.6166 

.6166 

3 

.6962 

.6963 

.6861 

.6861 

4 

• 7512 

.7513 

.7430 

.7430 

5 

• 7965 

• 7964 

.7896 

.7896 

6 

.8332 

.8333 

.8277 

.8277 

7 

.8635 

.8636 

.8589 

.8590 

8 

.8882 

.8882 

.8845 

.8845 

9 

.9085 

• 9079 

.9054 

.9048 

10 

.9250 

• 9258 

.9226 

• 9233 

ll 

.9388 

• 9396 

,9366 

• 9367 

12 

.9^97 

• 9497 

.9481 

• 9481 

13 

.9589 

.9589 

• 9575 

• 9575 

l4 

.9664 

.9664 

.9652 

• 9653 

15 

.9724 

.9724 

• 9715 

• 9715 

16 

• 9774 

• 9774 

.9767 

.9767 

17 

.9816 

.9815 

.9809 

.9809 

18 

.9848 

.9849 

.9844 

.9844 

19 

.9876 

.9876 

.9872 

.9872 

20 

.9899 

.9899 

• 9895 

• 9895 




TABLE II.- COMPARISON OE EXACT AND APPROXIMATE 
VALUES CO? TEE RESPONSE TO THE JUNCTION 


G-(s) = sin -2-s 
24 
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TABLE HI.- COMPARISON OP TEE RESULTS OP THE NUMERICAL CALCULATIONS 
OP THE FORCING- FUNCTIONS WITH THE EXACT VALUES 





fhponaeto arbitrary 

Arbitrary functionals) /ndbai response, Acs) forcing function, A(s) 
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S-z S/ O $/ ■$£ 

77/ne \ar/at/e, s 

Figure J . - Unit jump, grad/ent, and parabolic functions. 



<J/ $2 S3 Sq. 

v$ 


Figure . 2 . - Arbitrary function^ indicia!. response, and response 
to arbitrary function. 





Abstract 


Two matrix methods for calculating forcing 
functions from known responses are presented, one 
consisting of a numerical evaluation of Duhamel's 
integral and one consisting of straight-line or 
parabolic— arc approximations to the forcing function. 
The methods are suitable for application to many 
dynamic and some static problems. 
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